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Abstract 


We propose an efficient multistage method for solving a system of linear 
and nonlinear Volterra integral equations of the second kind. This numer- 
ical method is based on the Gauss—Legendre quadrature rule that obtains 
several values of unknown function at each step, and it will be shown that 
the order of the convergence is O(M~“), where M is the number of the 
nodes in the time discretization. The method has also the advantages of 
simplicity of application, less computational time, and useful performance 
for large intervals. In order to show the efficiency of the method, numerical 
results for the avian human influenza epidemic model is obtained that is 
comparable with the fourth-order Runge-Kutta method. 
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1 Introduction 


Numerical methods for solving a variety of system of Volterra integral equa- 
tions (SVIE) are offered in many publications. Here readers may refer to 
(1, 3, 7, 10] and others. 

In this article, we introduce a multistage numerical method for solving SVIE, 
and then the introduced method will be used to the avian human influenza 
epidemic model, which recently has been described by Iwami, Takeuchi, and 
Liu [6]. This mathematical model is proposed to interpret the spread of the 
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avian influenza from the bird world to the human world. The avian influenza 
has caused the death of millions of birds (almost 100% death) and it has 
become a disease of great importance both for animal and human health. 
Fortunately, there is still no evidence that the avian influenza virus can be 
transmitted among humans [12], but this virus is unstable and lacks of gnomic 
proofreading mechanism. Therefore the small errors, which occur when the 
copies themselves, go undetected and uncorrected. Since the specifics of mu- 
tations and evolution of the influenza viruses cannot be predicted, it is hardly 
possible to know if or when a virus such as the avian influenza, might ac- 
quire the properties needed to spread easily and sustainably among humans 
[5]. However, experts warn about an occurrence of the so-called “mutant 
the avian influenza”, which can be easily transmitted among humans with 
potentially devastating consequence. The US Congressional Budget Office 
has formally modeled the likely consequences of pandemic influenza and es- 
timates that up to 2 million of the US population might die, with up to 40% 
of all workers ill for as long as three or more weeks [11]. Thus, it is necessary 
and urgent to investigate the transmission process of the avian influenza and 
to take efficient measures to control the spread of the avian influenza. 


2 Preliminaries 


In this section, we introduce notations, definitions, and preliminaries facts, 
which are used throughout this article. Consider SVIE 


Vi) = G0) + [KU 5.Y)as On es22 7, (1) 


where Y(t) = [y™ (t), y?)(t),..., y (b]* is the desired function, 
G(t) = [9™ (#), 9 @),..., gM], and 


K(t,3,¥%(s))=|. 0 a ; 
kn(t, 8,y™ (t), y (0), ..-, 9) 


are the source functions (A‘ is the transposed matrix of A). Let system (1) 
is uniquely solvable. Necessary and sufficient conditions for the existence 
and uniqueness of solution (1) can be found in [9]. Therefore we assume the 
following conditions: 

(I) The source function G(t) is continuous (i.e., each component is continu- 
ous). 

(II) The kernel K(t,s,U) is a continuous function for 0 < s <t < T and 
0 <|| U [|< 0. 
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(III) The kernel K satisfies the Lipschitz condition 
| K(t,s,U) —K(T,s,V) |< L||U-V |, 


where the norm is defined as || U ||= max |u,(t)]. 
O<i<n 


3 Numerical procedure 


Let I, = {tm = mh,m = 0,1,..., MW} with ty, = Mh = T bea given uniform 
mesh on [0,7]. For given real numbers c; with 0 = cp < c1 < +--+ < cp =1, we 
choose the set of mesh points as Tp, := {tm,j = tm +cjh,m =0,1,...,M— 
1,7 =0,1,...,p}, where c;, 7 =1,2,...,p—1, are the roots of the (p — 1)th 
Legendre polynomial. The Legendre polynomials P,-1(x) are orthogonal 
with respect to the weight functions w(a) = 1 over the interval [—1, 1], that 
is, 


1 
/ Pral)Pa (de = Sens 
-1 


where Omn is the Kronecker delta (see, for example, [8, 4] for Legendre poly- 
nomials and [13] for funding roots of polynomials). 

For simplicity, let p = 3 (ec, = avs CQ= 3+V3) and let n = 2 (system 
dimensional). Then we will have a multistep method with at least four order 
of convergence, which one can generalize it by increasing the Legendre poly- 
nomial degree. 

Assume that y, & yO(tmj), m = 0,1,...,M—1, j = 0,1,2,3 and that 


Ying © Y (tm) = [y™ (tm,j), y (tm.z)]’- Thus system (1) can be rewrite as 
tm 

y (Ema) =I tng) +f” kaltmgs 8.0 (8)s¥(s))a 
0 


tm, j 
+ ka (tmajs 8 y(8), y(s))ds, (2) 
t 


tm 
ym§) = ms) +f ka(tmajs 8 yD (s), y(s))ds 


bm, j 
+f ka(tm,js 8, y) (s), yO s))ds, (3) 
t 


m 


m=0,1,...,M—1, j=1,2,3. 


Note that Yoo = dodo! = [g™ (0), g)(0)]', and suppose that approxi- 
mations Yo 1, Yo,2, Y1,0,---,Yn,o, have been calculated in the previous steps. 
Therefore the first integral in relations (2) and (3) can be estimated by using 
the Gauss—Legendre quadrature rule, 
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ki(tm,j,8,y" N(s),y (2) (s pie fs (ti,jy 8s y?(s),y(s))ds 


[ra(t m,j? tory), yO) + krlt bmg tof. ¥)| = Ai(m, j), 


[Fa(tm.js tins (2 af + ese =: A(m, j). 


(4) 


If 7 = 3, then the second integral in relations (2) and (3) can be similarly 
approximated as 


tm,3 
i: kr (tm,3,8,9(s), y(s))ds 
t 


m™m 


h 
mee tate alec) + iy (ten as tes Uns Yara) , 


tm,3 
/ ka(tm,3,8,9(s),y(s))ds 
t 


m 


h 
~ 2 [ks (tm,3, tm,1) yas y?),) + ke (teres tm,25 yas y)| ’ (5) 


and for 7 = 1, by using this quadrature rule, we obtain 


tm,1 
/ kr (tmas8;9(s), y(s))ds 
t 


m 


ch ha h 
a | tadcta tA Ga) Ga = 
5 ia Ly + ey ( + @)y ( = 5) 


2- 2- 2- 
ths tant + 2h, 9 (tm + 2h) yl + On| , 


2 


tm,1 
i, ko(tmas8,9(s),y(s))ds 
t 


cyh 


h h 
y k mM, tm = (1) m = (2) tm = 
a ate ttm + Ey ‘(tm + E) Yr (tm + &)) 


2 2— 2= 
“ho (brats ton - au y) (tm 7 2= V8 yyy (ta “i a) : 


b 


Finally, for 7 = 2, we obtain 
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tm,2 
/ ki (tm,2;8,y(s),y(s))ds 
tm 
ah 


h h h 
ae ce Wi ites, cel 
2 u(tm.astm + 674 (tm + @)o¥ (tm + &)) 


+k, (tm,2; tm + 


2 2 
2+ V5 yO (tm + ER) Mtn + 248) ) 


tm,2 
: ka (tm,28,9(s), y (s))ds 
t 


m 
coh 


2 


h 
5) 


2 2 
2+ VB YO (tm + ER) Mtn + 2) : 


h 
atm. tm “ Fu bm aT ae (tin + 


+kg (tm,2; tm + 


(7) 


In relations (6) and (7), tm + 24V3 p and tm + # do not belong to the mesh 


points II;,. Then we will have a problem to computing y (tm, + 24V3 p) and 


y (tm + By, i = 1,2. In order to overcome to this difficulty, we use the 
Lagrange interpolation 


for © =tm + 2£V3 ph or ¢ — tm + A Substituting these approximations into 
(2), leads to the system of algebraic equations 


; : ch h h h 
ra =g (tm.1) + Ai(m, 1) + a Astmastm oe G Piltm a G): Paltm A ra) 


2-—V3 2—V3 2-—V3 
= Pi (tin + “ Po (tm + oa 


+k; (tm,1, tm + 


h)) 


- 


j ‘ coh h h h 
Ys =9 (tm.2) + Ai(m, 2) + S* itm astm + 5, Paltm + 2); Paltm + 2) 


6 
9 44/3 9 44/3 2 aa 3 
+8 1, Putin + +8 5), Paltin + amt 


+kj (tm,2, tm + 


h)) 


? 


nee =9' (tm,3) + Ai(m, 3) 


h 
a D} [Fi (tm,3; tm,1s yor Wey) za Kilbenias tm,25 en ye)| ) 
po. (8) 
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This system can be solved by using an iterative method such as Newton— 
Raphson method or by using a suitable software package such as Maple or 
MATLAB 


4 Convergence analysis 


Theorem 1. The approximation method given by system (8) is convergent, 
and its order of convergence is at least 4 for the functions k; and g (i = 1, 2) 
with at least fourth-order continuous derivatives. 


Proof. Define Em,j = Yin,j — Y(tmj), where Y(tm,j),Ym,j € R? denote, 
respectively, the exact and approximate solutions of (1) at the point t = tm,;. 
Subtracting (6) from (2) and (3) (for 7 = 1) leads to 


Em =Yn1 — Y (ta) 
h m—l1 
=5 [K (tm,11ti,15 Ye) + K(tm,1, ti,25 Yio) 


h h 
+ F*|K(tmartm + 31 Plltm + €)) 
2-—v3 2-—VvV3 
81, Pltn £ iy] 


tis tm,1 
-f K (tm.1,8, Y*(s))ds -f Kia ONG. 
0 tin, 


+ K(tma;tm + 


where P(a) = [P1(x), Po(x)|. By adding and diminishing the terms 


BR 


m— 


h 
5 [K (tm,1,ti,1, ¥°(ti1)) + K (tm,1, ti,2, Y"(ti2))] , 
0 


+= 


and 


2-3 


h h h 2—V/3 
> tmnt + gry (tm + 6”) + K(tma,tm + LPG -+ 


we can write 


m—-1 


[Emil Sz Do [NK Ctm,rstias Yea) — K (tm, tia, ¥"(tia))| 
i=0 


+||.K (tm,15 4,2) Yi'n) — K (tm,1, ti,2, Y*(ti,2))|| 
h 


h h 
g)) ~ Kltmartm + EY (im + @ II 


cyh h 
a Kk tm stm Rp? tm 
+B [1K Utmastm + Bs Plém + : : 
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Di 4/3 2-4/3 
31, Pltm + as 


arts, ie oe “ 


+ ||K(tmjaytm + 


2 
—K(tma, tm 5 


+ (m+ Ille(@)II- 


In the previous relation, e(G) is the upper bound of the Gauss—Legendre 
integration error, that is, 


4 


ch h h 2—V3 2—V3 
elle (Knastm t ge Y (tm t )) + K(tm,i,tm + YS KY tn 1) 


fs m—-1 re 
NS So [mas tia, ¥*Ci)) + Kltmastia, ¥(t2))] -yf K(tm,a,8,¥"(s))ds|| 
i=0 i=0 % ti 


6 


tm, 
-| K (tm,1, 8, ¥*(s))ds|| < m|le(G)|| + |le(@) ||. 
ten 
If kj € C4((0, 1]) (¢ = 1,2), then e(G) = zagh*k (t,t, Y*(é)) where 0 <t < 
1; see [2]. 

In the following, by using the Lipschitz condition for the kernel K, we have 


Ball ES (all + Bal) + [Pen +) ¥en + PDI 
i=0 
+ ||Pltm += ny) iG <5) + (m+ I)lle(@)II- 


Let I(t) be the Lagrange interpolation error, that is, I(t) = Y‘(t) — P(t). 
Thus 


hL m-1 
lEmall < aa (|Zial] + || £i,2I1) 
1=0 
CAL h 2-3 
+ [|Z(tm + ZIM + |Z (tm 4 h))||} + (m + 1)Ile(G)|I. 
2 6 6 
Without loss of generality, assume max, \|Ems|| = ||/Em||. Then it is easy 
j=1,2, 
to see that 
m-1 
||Em,il] < AL SS ||Eial| + evhL||Z|| + (m+ Yle(@)II, 
i=0 


where J is the maximum error of the Lagrange interpolation. Finally, from 
the Gronwall inequality [2], we conclude that 
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||Em,a|| < (cA L||Z]| + M|le(G)||) exp(Ltm). 


This yields ||En,1|| + 0 as h > 0 and for the functions ky, k2, and g with 
at least fourth-order continuous derivatives, we have e(G) = O(h*) and I = 
O(h*). So we expect the error generally to be at least O(h*) as the numerical 
results confirm it. 


5 Application 


5.1 Description of the model 


In this part, the introduced numerical method is applied to interpret the avian 
human influenza epidemic model. Iwami, Takeuchi, and Liu [6] proposed the 
following mathematical model to interpret the spread of the avian human 
influenza epidemic 


X'=n- dX —wXY, 

Y’ =wXY — (b+ m)Y, 

S! =- pS — (AY - &H)S, 

B’ = Bi SY —(w+di +.)B, 

H! = 62.SH + «B-(w+a+y)H, 

R'=7H — pR. (9) 


All birds and humans in the effective population are divided into six main 
groups, respectively, including susceptible birds (X), birds infected with the 
avian influenza (Y), susceptible humans (S$), humans infected with the avian 
influenza (B), humans infected with mutant the avian influenza (H), and 
humans recovered from mutant the avian influenza (R). It is assumed that 
all birds infected with the avian influenza are dead or remain infected and 
can never recover. The parameter 7 and are the birth rates for birds and 
humans, respectively. The natural death rate is assumed to be b and yu for 
birds and humans, respectively. Furthermore, m and d are the death rates 
infected by wild the avian influenza, respectively, for birds and humans, and 
a is the additional death rate induced by mutant the avian influenza. The 
parameters w and /; are the rate at which the avian influenza is contracted 
from an average infected bird, (2 is the transmission rate of mutant the avian 
influenza in humans, 7y is the recovery rate, and ¢ is the mutation rate. 

In [6], it was investigated mathematical properties and qualitative anal- 
ysis of the above model, also was defined the basic reproductive number for 
infected birth with the avian influenza (“r9, which is the number of newly 
infected birds that are produced from anyone infected bird when all birds are 
susceptible”) and the basic reproductive number for infected humans with 
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mutant the avian influenza (Ro) as 


NW 


b(b+m)’ 


AB2 
wWetaty) 


Finally, using some theorems, global analysis to spread of the avian influenza 
and mutant influenza in the human world is obtained as follows: 

“If ro < 1 and Ro < 1, then the avian influenza and mutant the avian 
influenza cannot spread in the human world. On the other hand, if Ro > 1, 
then mutant the avian influenza spreads in the human world. Moreover, if 
ro > 1, then both the avian influenza and mutant the avian influenza spread 
in the human world.” 

It must be pointed that all constants are positive in the system (9) and that 
€ is sufficiently small. Moreover, b is sufficiently larger than pz (b > pz), a is 
less than d, and d is less than m (m > d > a) because of the differences of 
the virulence [6]. 


Ro = 


To = 


5.2 Solution procedure 


One can transfer the system (9) to the SVIE 
XOH= XO tre ‘i (b + w¥(s))X(s)ds, 
Y(t) = Y(0) +f (wX(s) —b—m)Y(s)ds, 
0 
S(t) = $(0) + — f (w+ ¥(s) + HH(s))S(s)ds, 


Bit) = BO) + | (G15(s)Y(s) —(u+d+ €)B(s))ds, 


H(t) = H(0) + | (B25(s)H(s) + <B(s) — (uta+7)H(s))ds, 
R(t) = R(0) + | (yH(s) — pR(s))ds. (10) 


In this way, we can apply the multistage method for this SVIE. In order to 
do it, set t= t,,,. Therefore 


X(t) =X(0) +n — f “(b+ w¥(s))X(s)ds— f (epee alae: 
Y (tn) =¥ (0) + i: " (wX (8) —b a) ¥ (ade 4 | (a —b rede, 


n 
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46 
tn 
S(tp,4) =5(0) + Ming — i: (u + AY (s) + BoH(s))S(s)ds 
0 
tni 
= i! (u + BLY (s) + B2H(s))S(s)ds, 
tn 
tn 
Bltns) =BO)+ ["(B.S()¥(s) ~ (w+ d+ QB(s))as 
0 
tnyi 
+ [ @S()¥() ~ ut d+ )B(s))ds, 
th 
tn 
A(tys) =H(0) + | (895(s)H(s) + €B(s) —(u+aty)H(s))ds 
0 
thi 
+f (825(s)H(s) + €B(s) — (u+taty)H(s))ds, (11) 
tn 
tiny tri 
Rtas) =RO0)+ f" OH(s)— wR(s))as + f™ OH(s) ~ HR(s))as 
G20 G ca <4 HAS. 
Assume that Xn, Ynji,---, Rn, are approximation for exact solutions X (¢), 
Y(t),..., R(t) in the point t,,;. Similar to the previous section, we can ap- 


proximate the integrals in the SIVE (11). Use two points Gauss—Legendre 
for integration on the interval [0,t,,]. Thus 


é - 
he 

[ (b+ wY(s)) swt [((b + w¥51)X5,1 + (0+ wYj,2) X52] =: Ar, 

0 j=0 


n-1 
h 
mS DD [(wXj,1 —b-— m)Yj,1 + (WX5,2 —b-— m)Yj,2] =: Aa, 
2m 
ts h n-1 
/ (yH(s) — wR(s))ds » 5 Yo yy — Ry + 1A, — wR32] =: As- 
0 j=0 


Continuing the same process for other integrals leads to the system of alge- 
braic equations 


ev | 7 ee | 


hey 
h))Pi(tn + 


2 


(b+ wPo(tr + h) 


Xn =Xo+ nin — At 


+(b+ wPa(tn + = 7) )Palte +3), 
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Xn2=Xo+Ntn2—- At 


— 92 (b+ wPalta +2 +S A yPilty re oy 
+(b + WP o(tn + *))PuUtn + *) 5 
Xn3 =Xo + Htn3 — At 
om [(b+ wn )Xna + (b+ wYn2)Xn2], 
Yn, =Yo + Ag 
+ no (wPiltn + = 8 b= m)Paltn += ae) 
+(wP1 (tr + ) —b-—m)Pi (tn + al , 
Yn,2 =Yo + Ag 
+ nea (WP1 (tn + : Sy — b= m)P2(tn + z ual 
+UPAltn + §)— b= m)Pulln + 5), 
Yn,3 =Yo + Aa 
- [(wXna —b- m)Yna + (WXn,2)¥a,al 
Rn =Ro + Ao 
+ no *Psltn + Si) — [UP oltn + =i) 
+1Poltn +5) — HPalln +B) 
Rn2 =Ro + As 
+72 | oPs(tn +7 amc See aA) 


h h 
+1Poltn + 5) ~HPaltn + 5) 
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Rn3 =Ro + Ab 
h 
~ 5 Hn — HRna + Hn — HRn,2], n=0,1,...,.N—1, 
where 


3 
X (tn tha) © Pi(tn + ha) = S° L; (tr +h2)Xn,j, 


Y (ty + ha) ye Pa(tn are ha) = L; (tn 9 hz)Yn,j; 


3 
R(ty + ha) © Pe(tn + hax) = YG in) Rigg, 
j=0 


3 
for the Lagrange polynomial L(t, + hx) = [[ = and gz = 2e/3 or 
m=0 
mA; 
r= 1, 
In each step (n = 0,1,..., N —1), three values of unknowns functions can be 


obtained by solving this system. 


6 Numerical results and discussion 


In this section, we present some numerical results to investigate the spread 
of the avian influenza. The initial value of system (10) is fixed at X(0) = 10, 
Y(0) = 2, S(0) = 100, B(O) = 0, H(0) = 0, and R(0) = 0. These amounts 
are chosen from the region 2 := {(X,Y,5,B,H,R): xX >0,Y >0,S > 
0,B =0,H = 0, R = 0}, which denotes that there do not exist the infected 
humans with the avian influenza and mutant the avian influenza. 

Similar to [6], the parameters are fixed at b = 5, m = 5, A = 3, wp = 0.015, 
B, =0.2,€ = 10-3, d= 1, a = 0.06, and y = 0.01. We will provide numerical 
results of infected birds (Y) with the avian influenza, infected humans (B) 
with the avian influenza, and infected humans (H) with mutant the avian 
influenza by choosing different constants 7, w, and fo. 

The numerical experiments are carried out in Maple, and nonlinear systems 
arising from the nonlinear integral equations are solved by Maple routine 
fsolve. 


Example 1. In system (10), let 7 = 26.5, w = 2 and 82 = 0.003. Then 
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nwo Bor 


— == Sie, Bee wr 
b(b +m) * peat) 


To = 


Figures 1 and 2:(a) describe the pandemic effect on the birds and the 
human world. Figure 1:(a) describes an endemic in the infected birds just 
after the occurrence of the avian influenza. According to the reported global 
analysis from [6], Figure 1:(b) shows an initially pandemic of the avian in- 
fluenza and afterward pandemic level would be decreased, and Figure 2:(a) 
interprets that the infected humans with mutant the avian influenza outbreak 
and eventually keep the relatively high level of the size. It suggests that the 
mutant the avian influenza pandemic will occur if we do not take efficient 
measures to control the spread of the avian influenza. Moreover, CPU times 
of the computations are shown in the figures. 


(a) (b) 


25 


05 


Figure 1: (a) The bird rate infected with avian, T = 10, M = 100, CPU time=19.5s. (b) 
The human rate infected with the avian influenza, T = 100, M = 1000, CPU time=408s. 
ro = 1.06, Ro & 7.1. 


Example 2. The parameters are the same as the previous example except for 
(G2 = 0.0015 (the transmission rate of mutant the avian influenza is reduced 
by half). Therefore the basic reproduction number for the avian influenza 
in the bird and human world, respectively, is equal to rg = Weim) = 1.06, 


— fer 
Ro = nutary © 3.5. 


Figures 2:(b) and 3:(a), respectively, show the endemic process in the bird 
world after the morbidity of the avian influenza and the occurrence of the 
pandemic of the avian influenza in the human world. Figure 3:(b) shows the 
human rate infected with mutant the avian influenza. Since Ro is relatively 
small, we do not have mutant the avian influenza pandemic and in this case, 
we have only the avian influenza pandemic. 
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HO Yo 


a 


Figure 2: (a) The rate infected with mutant the avian influenza, T = 100, M = 1000, 
CPU time=408s, ro = 1.06, Ro © 7.1. (b) The bird rate infected with the avian influenza, 
T = 10, M = 100, CPU time=19.3s, ro = 1.06, Ro © 3.5. 


Example 3. Let 7 = 30 and w = 1.5; then rg = 0.9 and Ro & 7.1. 


The numerical results for this example are shown in Figures 4 and 5. In 
these figures, there are two pandemic for the avian influenza and mutant the 
avian influenza even, when all the infected birds and humans with the avian 
influenza were extinct. 


(a) (b) 


Figure 3: (a) The human rate infected with the avian influenza, (b) The rate infected 
with mutant the avian influenza. T = 100, M = 1000, CPU time=398s. ro = 1.06, Ro © 
3.5. 
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(a) 


Yo 


Os 


Figure 4: (a) The bird rate infected with the avian influenza, T = 10, M = 100, CPU 
time=18.2s. (b) The human rate infected with the avian influenza, T = 100, M = 1000, 
CPU time=409s. ro = 0.9, Ro & 7.1. 


From these numerical results similar to [6], we conclude the following 
result: 


e Elimination policy of the infected birds is not necessarily useful, because 
the speed of spread of disease in the bird world is fast and their lifetime 
is shorter than humans. 


e For ro < 1, the infected humans extinct, but there are some infected 
humans for ro > 1 (see Figures 1:(a), 3:(a), 4:(b)). 


e If Ro is not large enough, then the outbreak of mutant influenza will 
not occur (see Figures 2:(a), 3:(b), 5). 


The numerical results for Examples 1-3 are comparable with the fourth-order 
Runge-Kutta method [6] in the whole intervals. Moreover, the presented 
method, in contrast to Runge-Kutta methods, did not need to computation 
a lot of derivations, and the computation time for this method is short. Of 
course, we can increase the convergence order of the method. Also, numerical 
results are obtained for each time interval that is one of the other advantages 
of the introduced method. 
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Figure 5: The rate infected with mutant the avian influenza, T = 100, M = 1000, CPU 
time=409s, ro = 0.9, Ro & 7.1. 
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